{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"right\" style=\"max-width: 200px; height: auto\" src=\"cfds_logo.png\">\n",
    "\n",
    "###  Lab 07 - \"Unsupervised Machine Learning k-Means Clustering and Expectation Maximization\"\n",
    "\n",
    "Chartered Financial Data Scientist (CFDS), Spring Term 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last lab you learned about how to utilize **supervised** learning classification techniques namely (1) the Gaussian Naive-Bayes (Gaussian NB) classifier, (2) the k Nearest-Neighbor (kNN) classifier and (3) the Support Vector Machine (SVM) classifer. \n",
    "\n",
    "In this lab we will learn about two **unsupervised** machine learning technique referred to as **k-Means Clustering** and **Expectation-Maximization (EM) Clustering**. We will use this technique to classify un-labelled data (i.e., data without defined categories or groups). In general, clustering-based techniques are widely used in **unsupervised machine learning**.\n",
    "\n",
    "<img align=\"center\" style=\"max-width: 500px\" src=\"machinelearning.png\">\n",
    "\n",
    "The **k-Means Clustering** algorithm is one of the most popular clustering algorithms used in machine learning. The goal of k-Means Clustering is to find groups (clusters) in a given dataset. It can be used (1) to **confirm business assumptions** about what types of groups exist or (2) to **identify unknown groups** in complex data sets. Some examples of business-related use cases are:\n",
    "\n",
    ">- Segment customers by purchase history;\n",
    ">- Segment users by activities on an application or a website;\n",
    ">- Group inventory by sales activity; or,\n",
    ">- Group inventory by manufacturing metrics.\n",
    "\n",
    "(Source: https://www.datascience.com/blog/k-means-clustering)\n",
    "\n",
    "Once the algorithm has been run and the groups are defined, any new data can be easily assigned to the correct group."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, pls. don't hesitate to ask all your questions either during the lab or send us an email (using our\n",
    "fds.ai email addresses)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lab Objectives:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After today's lab, you should be able to:\n",
    "\n",
    "> 1. Understand how a **k-Means Clustering** algorithm can be trained and evaluated.\n",
    "> 2. Know how to select an **optimal number of clusters** or cluster means.\n",
    "> 3. Understand how an **Expectation-Maximisation (EM) Clustering** algorithm can be trained and evaluated.\n",
    "> 4. Know how to Python's **sklearn library** to train and evaluate arbitrary classifiers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup of the Jupyter Notebook Environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to the previous labs, we need to import a couple of Python libraries that allow for data analysis and data visualization. We will be using pandas, numpy, sklearn, matplotlib and seaborn library to conduct this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the pandas data science library\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# import the scipy spatial distance capability\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "# import sklearn data and data pre-processing libraries\n",
    "from sklearn import datasets\n",
    "from sklearn.datasets.samples_generator import make_blobs\n",
    "\n",
    "# import sklearn k-means and gaussian-mixture classifier library\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "# import matplotlib data visualization library\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "# import matplotlibs 3D plotting capabilities\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create nice looking plots using the **seaborn** plotting theme:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('seaborn')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enable inline Jupyter notebook plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppress potential warnings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set random seed of all our experiments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_seed = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. k-Means Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Dataset Download and Data Assessment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The iris dataset is a classic and very easy multi-class classification dataset. This dataset consists of 3 different types of irises’ (classes),  namely Setosa, Versicolour, and Virginica) as well as their respective petal and sepal length (features). In total, the dataset consists of **150 samples** (50 samples per class) as well as their corresponding **4 different measurements** taken for each sample. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"center\" style=\"max-width: 700px; height: auto\" src=\"iris_dataset.png\">\n",
    "\n",
    "(Source: http://www.lac.inpe.br/~rafael.santos/Docs/R/CAP394/WholeStory-Iris.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please, find below the list of the individual measurements (features):\n",
    "\n",
    ">- `Sepal length (cm)`\n",
    ">- `Sepal width (cm)`\n",
    ">- `Petal length (cm)`\n",
    ">- `Petal width (cm)`\n",
    "\n",
    "Further details on the dataset can be obtained from the following publication: *Fisher, R.A. \"The use of multiple measurements in taxonomic problems\" Annual Eugenics, 7, Part II, 179-188 (1936); also in \"Contributions to Mathematical Statistics\" (John Wiley, NY, 1950).\"*\n",
    "\n",
    "Let's load the dataset and conduct a preliminary data assessment: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = datasets.load_iris()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print and inspect feature names of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.feature_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print and inspect the class names of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.target_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print and inspect the top 5 feature rows of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(iris.data).head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print and inspect the top 5 labels of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(iris.target).head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine and print the feature dimensionality of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine and print the label dimensionality of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris.target.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's briefly envision how the feature data is collected and recorded in the data:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"center\" style=\"max-width: 900px; height: auto\" src=\"featurecollection.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data distributions of the distinct features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 10))\n",
    "iris_plot = sns.load_dataset(\"iris\")\n",
    "\n",
    "# supervised scenario\n",
    "# sns.pairplot(iris_plot, diag_kind='hist', hue='species');\n",
    "\n",
    "# unsupervised scenario\n",
    "sns.pairplot(iris_plot, diag_kind='hist');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now imagine that we did not have the convenient `species` (class) label associated with each observation in the iris dataset. **How could we distinguish the three species in the dataset?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. The k-Means Clustering Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As already mentioned the **k-Means Clustering** algorithm is one of the most popular \"first choice\" unsupervised clustering algorithms to find groups (clusters) in a given multidimensional dataset $X$.\n",
    "\n",
    "<img align=\"center\" style=\"max-width: 400px; height: auto\" src=\"kmeans.png\">\n",
    "\n",
    "Thereby, the basic form of k-Means Clustering makes the following **two assumptions**:\n",
    "\n",
    "- Each observation is closer to its own cluster center than to the center of the other clusters.\n",
    "- A cluster center is the arithmetic mean of all the points that belong to the cluster.\n",
    "\n",
    "Let's briefly revisit the distinct step of the algorithm before applying it to the iris dataset. Therefore, let's assume:\n",
    "\n",
    "- We have dataset $X$ consisting records $x_1, x_2, x_3, ..., x_n \\in \\mathcal{R}^d$; \n",
    "- That samples are clustered around $k$ centers (the \"$k$ means\") denoted by $\\mu_1, \\mu_{2}, ..., \\mu_{k} \\in \\mathcal{R}^d$; and,\n",
    "- Each sample $x_{i}$ belongs to its closest mean $\\mu_{i}$.\n",
    "\n",
    "We can then iteratively perform the following steps that comprise the **k-Means Clustering** algorithm:\n",
    "\n",
    ">- **Step 1** - Pick $k$ random points $\\mu_{i}$ as cluster centres called 'means'.\n",
    ">- **Step 2** - Assign each $x_i$ to its to nearest cluster mean by calculating its distance to each mean.\n",
    ">- **Step 3** - Determine the new cluster centres by calculating the average of the assigned points in each cluster.\n",
    ">- **Step 4** - Repeat Step 2 and 3 until none of the cluster assignments change.\n",
    "\n",
    "Note, that a single execution of all the four steps outlined above is usually referred to as 'iteration'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3. k-Means Clustering in a 2-Dimensional Feature Space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's see how we can apply it to the iris dataset. We will start with an introductory example of detecting the classes of the iris dataset based on two of its features namely the (1) `Petal length (cm)` and (2) `Petal width (cm)`. Let's first gain an intuition of those two features as well as their distribution by visualizing them accordingly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length (3rd feature in the dataset) vs. petal width (4th feature in the dataset)\n",
    "ax.scatter(iris.data[:,2], iris.data[:,3])\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[petal_length]\", fontsize=14)\n",
    "ax.set_ylabel(\"[petal_width]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Petal Length vs. Petal Width Feature Distribution', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now define the parameters of the k-Means Clustering. We will start by specifying the **number of clusters** $k$ we aim to detect in the iris dataset. We hypothesize that our observations are drawn from an unknown distributions of three iris flower species (each distribution corresponding to a different mean $\\mu_1$, $\\mu_2$, and, $\\mu_3$). Therefore, we set the number of clusters to be detected to $k=3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_clusters = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to define a corresponding number of **initial 'means' $\\mu_{i}$** (the initial random cluster centers) that will be used as 'starting points' in the first iteration of the clustering process. In our case we will specify $k=3$ cluster means each of dimension 2, since we aim to retrieve 3 clusters based on the 2 features `Petal length (cm)` and `Petal width (cm)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_means = np.array([[1.0, 3.0], [2.0, 6.0], [1.0, 7.0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will define a **maximum number of iterations** that we want to run the k-Means Clustering algorithm. Please, note that the clustering terminates once there will be no further changes in the cluster assignments. However, it's good practice to define an upper bound of the iterations applied in the clustering (especially when analyzing datasets that exhibt a high-dimensional feature space):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_iterations = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we ready to initialize an instance of the **k-Means Clustering** algorithm using Python's `sklearn` library of data science algorithms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=no_clusters, init=init_means, max_iter=max_iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the k-Means Clustering to learn a model of the `Petal length (cm)` and  `Petal width (cm)` features. Pls. note that we are using columns 2 and 3 to extract the values of the two features from the iris dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans.fit(iris.data[:,2:4]) # note that we are using column 2 (petal length) and 3 (petal width) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have conducted the clustering, let's inspect the distinct cluster labels that have been assigned to the individual records of the iris dataset. This can be achieved by calling the `labels_` attribute of the fitted model: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = kmeans.labels_ # obtain the assigned cluster labels\n",
    "print(labels)           # print the cluster labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, we want to inspect the coordinates of the cluster means (sometimes also referred to as \"centroids\") assigned by the algorithm. This can be achieved by calling the `cluster_centers_`attribute of the fitted model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = kmeans.cluster_centers_ # obtain the assigned cluster means \n",
    "print(means)                    # print the cluster center coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now plot the iris dataset records using the two features `Petal length (cm)` and `Petal width (cm)` as well as the labels and cluster means determined by the **k-Means Clustering** algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length vs. petal width and corresponding classes\n",
    "scatter = ax.scatter(iris.data[:,2], iris.data[:,3], c=labels.astype(np.float), cmap=plt.cm.Set1)\n",
    "\n",
    "# prepare data legend\n",
    "legend = ax.legend(*scatter.legend_elements(), loc='upper left', title='Cluster')\n",
    "\n",
    "# add legend to plot\n",
    "ax.add_artist(legend)\n",
    "\n",
    "# plot cluster means\n",
    "ax.scatter(means[:,0], means[:,1], marker='x', c='black', s=100)\n",
    "\n",
    "# iterate over distinct cluster means\n",
    "for i, mean in enumerate(means):\n",
    "    \n",
    "    # determine max cluster point distance\n",
    "    cluster_radi = cdist(iris.data[:, 2:4][labels==i], [mean]).max()\n",
    "    \n",
    "    # plot cluster size\n",
    "    ax.add_patch(plt.Circle(mean, cluster_radi, fc='darkgrey', edgecolor='slategrey', lw=1, alpha=0.1, zorder=1))\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[petal_length]\", fontsize=14)\n",
    "ax.set_ylabel(\"[petal_width]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Petal Length vs. Petal Width - Clustering Results', fontsize=14);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%HTML\n",
    "<div align=\"middle\">\n",
    "<video width=\"60%\" controls>\n",
    "<source src=\"kmeansvideo.mp4\" type=\"video/mp4\">\n",
    "</video></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be observed that **k-Means Clustering algorithm** nicely found three clusters in the dateset. Let's inspect to which extend this corresponds to the true 'species' class labels 'verginica', 'setosa', and 'versicolor' of the iris dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 8))\n",
    "\n",
    "#### plot true iris class labels\n",
    "\n",
    "# add grid\n",
    "ax[0].grid(linestyle='dotted')\n",
    "\n",
    "# iterate over distinct species\n",
    "for species in np.unique(iris.target):\n",
    "    \n",
    "    # obtain iris petal length and petal width\n",
    "    iris_features = iris.data[iris.target == species,:]\n",
    "    \n",
    "    # obtain iris species name\n",
    "    iris_target_name = iris.target_names[species]\n",
    "    \n",
    "    # plot petal length vs. petal width as well as the true labels\n",
    "    ax[0].scatter(iris_features[:,2], iris_features[:,3], c='C{}'.format(str(species)), label=iris_target_name)\n",
    "\n",
    "# prepare data legend\n",
    "ax[0].legend(loc='upper left', title='Classes')\n",
    "\n",
    "# add axis legends\n",
    "ax[0].set_xlabel(\"[petal_length]\", fontsize=14)\n",
    "ax[0].set_ylabel(\"[petal_width]\", fontsize=14)\n",
    "\n",
    "#### plot clustering results\n",
    "\n",
    "# add grid\n",
    "ax[1].grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length vs. petal width and corresponding classes\n",
    "scatter = ax[1].scatter(iris.data[:,2], iris.data[:,3], c=labels.astype(np.float), cmap=plt.cm.Set1)\n",
    "\n",
    "# prepare data legend\n",
    "ax[1].legend(*scatter.legend_elements(), loc='upper left', title='Cluster')\n",
    "\n",
    "# plot cluster means\n",
    "ax[1].scatter(means[:,0], means[:,1], marker='x', c='black', s=100)\n",
    "\n",
    "# iterate over distinct cluster means\n",
    "for i, mean in enumerate(means):\n",
    "    \n",
    "    # determine max cluster point distance\n",
    "    cluster_radi = cdist(iris.data[:, 2:4][labels==i], [mean]).max()\n",
    "    \n",
    "    # plot cluster size\n",
    "    ax[1].add_patch(plt.Circle(mean, cluster_radi, fc='darkgrey', edgecolor='slategrey', lw=1, alpha=0.1, zorder=1))\n",
    "\n",
    "# add axis legends\n",
    "ax[1].set_xlabel(\"[petal_length]\", fontsize=14)\n",
    "ax[1].set_ylabel(\"[petal_width]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.suptitle('Petal Length vs. Petal Width - True Class Labels (left) and Clustering Results (right)', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, it seems that our clustering did a quite good job."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, let's inspect the distance of all dataset records $X$ to their nearest means $\\mu_{i}$. Let's recall that k-Means Clustering conducts a local optimization of the sum of \"squared errors\", as expressed by:\n",
    "\n",
    "$$E(\\mu_{1}, \\mu_{2}, ..., \\mu_{k}) = \\sum_{i=1}^{n}(x_{i}-\\mu_{k(i)})^{2},$$\n",
    "\n",
    "were $x_{i}$ denotes a feature vector (or observation) of the dataset and $\\mu_{k(i)}$ its closest mean in the feature space $\\mathcal{R}^{d}$.\n",
    "\n",
    "We can obtain the sum of those squared distances $E(\\mu_{1}, \\mu_{2}, ..., \\mu_{k})$ by calling the `inertia_` attribute of the learned k-Means Clustering model. It will return the sum of squared distances of each sample to its closest cluster center:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = kmeans.inertia_\n",
    "print(distances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We recommend you try the following exercises as part of the lab:\n",
    "\n",
    "**1. Train and evaluate the k-Means Clustering algorithm and obtain its assigned labels as well as squared distances $E(\\mu_{1}, \\mu_{2}, ..., \\mu_{k})$ for distinct max iterations.**\n",
    "\n",
    "> Continuously increase the number of training iterations $i$ of the k-Means Clustering starting with 1 and up to 5 iterations ($i=1,...,5$) and repeat the clustering accordingly. What can be observed in terms of the cluster means as well as the sum of squared cluster distances with increasing $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ***************************************************\n",
    "# INSERT YOUR CODE HERE\n",
    "# ***************************************************"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2. Determine if the k-Means Clustering algorithm always converges to the same result.**\n",
    "\n",
    "> Carefully review the k-Means algorithm and answer to the following question: Does the k-Means algorithm always converge to the same result? Please, explain your reasoning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ***************************************************\n",
    "# INSERT YOUR CODE HERE\n",
    "# ***************************************************"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3. Application of the k-Means Clustering algorithm to distinct data distributions.**\n",
    "\n",
    "> Consider the following data distributions. Determine which are suitable for a k-Means clustering and what $k$ value should be applied in the clustering. Please, explain your reasoning.\n",
    "\n",
    "<img align=\"center\" style=\"max-width: 600px; height: auto\" src=\"clustering.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ***************************************************\n",
    "# INSERT YOUR CODE HERE\n",
    "# ***************************************************"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4. k-Means Clustering in a 3-Dimensional Feature Space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enable interactive 3-dimensional plotting inside of Jupyter notebooks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from importlib import reload\n",
    "reload(plt)\n",
    "\n",
    "# import the seaborn plotting library\n",
    "import seaborn as sns\n",
    "\n",
    "# import matplotlibs 3D plotting capabilities\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now extend the k-Means Clustering to a 3-dimensional features space $\\mathcal{R}^{3}$ space by clustering the first three features, namely `Sepal length (cm)`, `Sepal width (cm)`, and `Petal length (cm)` of the Iris dataset. Let's again start by visually inspecting the 3-dimensional feature space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# init 3D plotting\n",
    "ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=30, azim=120)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length (3rd feature in the dataset) vs. petal width (4th feature in the dataset)\n",
    "ax.scatter(iris.data[:,0], iris.data[:,1], iris.data[:,2], s=40)\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[sepal_length]\", fontsize=14)\n",
    "ax.set_ylabel(\"[sepal_width]\", fontsize=14)\n",
    "ax.set_zlabel(\"[petal_length]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Sepal Length vs. Sepal Width vs. Petal Length', fontsize=14)\n",
    "\n",
    "# show the 3-dimensional plot\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start again by defining a **max. number of iterations** we aim to run the k-Means Clustering algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_iterations = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also initialize a corresponding number $k$ **initial random cluster 'means' $\\mu_{i}$**. This time will specify $k=3$ cluster means each of dimension 3, since we aim to retrieve 3 clusters based on the 3 features `Sepal length (cm)`, `Sepal width (cm)`, and `Petal length (cm)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_means = np.array([[1.0, 3.0, 3.0], [2.0, 6.0, 5.0], [1.0, 7.0, 2.0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we ready to initialize an instance of the **k-Means Clustering** algorithm using Python's `sklearn` library of data science algorithms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=no_clusters, init=init_means, max_iter=max_iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the k-Means Clustering to now learn a model of the `Sepal length (cm)`, `Sepal width (cm)`, and `Petal length (cm)` features. Pls. note that we are using columns 0, 1, and 2 to extract the values of the three features from the iris dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans.fit(iris.data[:,0:3]) # note that we are using column 1 (sepal length), 2 (sepal width) and 3 (petal length) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's again inspect the labels assigned to each individual record:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = kmeans.labels_ # obtain the assigned cluster labels\n",
    "print(labels)           # print the cluster labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as the determined cluster means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = kmeans.cluster_centers_ # obtain the assigned cluster means \n",
    "print(means)                    # print the cluster center coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also plot the iris dataset records using the three features `Sepal length (cm)`, `Sepal width (cm)` and `Petal length (cm)`, their corresponding learned labels as well as the cluster means as determined by the **k-Means Clustering** algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# init 3D plotting\n",
    "ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=30, azim=120)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length vs. petal width and corresponding classes\n",
    "ax.scatter(iris.data[:,0], iris.data[:,1], iris.data[:,2], c=labels.astype(np.float), cmap=plt.cm.Set1, s=40)\n",
    "\n",
    "# plot cluster means\n",
    "ax.scatter(means[:,0], means[:,1], means[:,2], marker='x', c='black', s=100)\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[sepal_length]\", fontsize=14)\n",
    "ax.set_ylabel(\"[sepal_width]\", fontsize=14)\n",
    "ax.set_zlabel(\"[petal_length]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Sepal Length vs. Sepal Width vs. Petal Length', fontsize=14);\n",
    "\n",
    "# show the 3-dimensional plot\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's again inspect the distance of all dataset records $X$ to their nearest means $\\mu_{i}$ by calling the `inertia_` attribute of the learned k-Means Clustering model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = kmeans.inertia_\n",
    "print(distances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enable inline Jupyter notebook plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5. Drawbacks of the k-Means Clustering Algorithm "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though the **k-Means Clustering** algorithm is one of the most popular clustering algorithms used in machine learning. There are three major drawbacks associated with it, which are explained hereafter:\n",
    "\n",
    "1. The k-Means Clustering is guaranteed to improve the result in each iteration but there are **no guarantees** that it will find the **global best of clusters**. Recall the example, where we initalize the algorithm with a different seed of cluster means. \n",
    "\n",
    "> **Practical solution:** Run the algorithm with multiple random initializations. This also done per default when using the `scikit` learn of machine learning algorithms.\n",
    "\n",
    "2. The k-Means Clustering **cannot learn** the **optimal number of clusters** from the provided data. E.g., if we ask the algorithm for 20 clusters it will find 20 clusters, which may or may not be meaningful. \n",
    "\n",
    "> **Practical solution:** Use the \"Elbow\" technique as explained in the next section of the notebook. Another option might be the usage of a more complex clustering algorithm like Gaussian Mixture Models, or one that can choose a suitable number of clusters, e.g., the DBSCAN clustering algorithm.\n",
    "\n",
    "3. The k-Means Clustering **doesn't work well** in instances of a **non-linear seperable** dataset. This is caused by its assumption that points will be closer to their own cluster center than to others.\n",
    "\n",
    "> **Practical solution:** Transform (if possible) the dataset into a higher dimension where a linear separation becomes possible, e.g., by using a spectral clustering algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.6. Optimal Cluster Number Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that, one of the basic ideas behind unsupervised machine learning methods, such as **k-Means clustering**, is to define clusters for which the total intra-cluster variation (usually measured by the total sum of squared distances) is minimized:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$k^{*} =\\underset{k}{\\arg \\min} \\sum_{i=1}^{n}(x_{i}-\\mu_{k(i)})^{2},$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "were $x_{i}$ denotes a single feature vector (or observation) in the dataset and $\\mu_{k(i)}$ its closest mean in the feature space $\\mathcal{R}^{d}$. Challenge: What is the optimal number of clusters $k$ for a given dataset? Selection of the right $k$ may result in the following issues:\n",
    "\n",
    "- if $k$ too small (under-segmentation), then the clusters are too diverse; and;\n",
    "- if $k$ too high (over-segmentation), then the clusters are too fine-grain.\n",
    "\n",
    "Examples: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"center\" style=\"max-width: 800px; height: auto\" src=\"kselection.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution: We can then use the sum of \"squared errors\" $E(\\mu_{1}, \\mu_{2}, ..., \\mu_{k})$ metric to find an optimal number of clusters $k$! This can be achieved by the execution of the so-called **'Elbow'** technique defined by the following algorithm:\n",
    "\n",
    ">- **Step 1** - Compute the k-Means clustering algorithm for different number of clusters $k$.\n",
    ">- **Step 2** - For each $k$ calculate the sum of the within-cluster sum of squared distances $E(\\mu_{1}, \\mu_{2}, ..., \\mu_{k})$.\n",
    ">- **Step 3** - For each $k$ plot the $k$ value vs. its corresponding sum of within-cluster sum of squared distances $E$. \n",
    ">- **Step 4** - Inspect the plot and determine the location of a bend (appropriate number of clusters)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's utilize the **'Elbow'** technique by first defining a max. number of iterations that we aim to apply at each k-Means clustering run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_iterations = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can implement the and run the 'elbow' technique:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the list of squared distances\n",
    "sum_of_squared_distances = []\n",
    "\n",
    "# define the range of k-values to investigate\n",
    "K = range(1,30)\n",
    "\n",
    "# iterate over all k-values\n",
    "for k in K:\n",
    "    \n",
    "    # init the k-Means clustering algorithm of the current k-value\n",
    "    kmeans = KMeans(n_clusters=k, init='random', max_iter=max_iterations)\n",
    "    \n",
    "    # run the k-Means clustering of sepal-length and sepal-width features\n",
    "    kmeans = kmeans.fit(iris.data[:,0:2])\n",
    "    \n",
    "    # collect the sum of within squared distances of the current k-value\n",
    "    sum_of_squared_distances.append(kmeans.inertia_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Upon completion of the loop above let's inspect the distinct within-cluster sum of squared distances $E$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print the collected sum of squared distances of each k\n",
    "sum_of_squared_distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, let's plot the cluster number $k$ vs. the within-cluster sum of squared distances $E$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length (3rd feature in the dataset) vs. petal width (4th feature in the dataset)\n",
    "ax.plot(K, sum_of_squared_distances)\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[number of clusters $k$]\", fontsize=14)\n",
    "ax.set_ylabel(\"[within-cluster distance $E$]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Cluster Number $k$ vs. Within-Cluster Distance $E$', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We recommend you to try the following exercises as part of the lab:\n",
    "\n",
    "**1. Apply the k-Means Clustering algorithm to all four features contained in the Iris dataset.**\n",
    "\n",
    "> Use the k-Means classifier to learn a model of all four features contained in the Iris dataset, namely `Sepal length (cm)`, `Sepal width (cm)`, `Petal length (cm)` and `Petal width (cm)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ***************************************************\n",
    "# INSERT YOUR CODE HERE\n",
    "# ***************************************************"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2. Determine the optimal number of cluster values $k$ of all four features contained in the iris dataset.**\n",
    "\n",
    "> Determine the optimal number of clusters $k$ needed to cluster the observations of all four features contained in the iris dataset using the **'Elbow'** technique outlined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ***************************************************\n",
    "# INSERT YOUR CODE HERE\n",
    "# ***************************************************"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Expectation Maximization (EM) Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"One way to think about the k-means model is that it places a circle (or, in higher dimensions, a hyper-sphere) at the center of each cluster, with a radius defined by the most distant point in the cluster\" (Jake VanderPlas, The Python Data Science Handbook). To investigate this further let's have a closer look at another very skewed data distribution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Dataset Creation and Data Assessment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a very skewed distribution that exhibits a eliptical characteristic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a sample data distribution\n",
    "data, labels = make_blobs(n_samples = 1000, centers = 3, cluster_std = .7, random_state = 0)\n",
    "\n",
    "# init the random state and skew the data\n",
    "random_state = np.random.RandomState(8)\n",
    "data_skewed = np.dot(data, random_state.randn(2, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once created, let's visualize the skewed distribution accordingly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot the two dimensions of the skewed distribution\n",
    "scatter = ax.scatter(data_skewed[:,0], data_skewed[:,1])\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[feature $x_1$]\", fontsize=14)\n",
    "ax.set_ylabel(\"[feature $x_2$]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Sample Skewed Distribution', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, the distribution looks indeed very skewed. It furthermore seem to consist of three commingled eliptical clusters. Let's see if we can identify those clusters using the k-Means Clustering algorithm that we discussed before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the k-Means Clustering algorithm\n",
    "kmeans = KMeans(n_clusters=3, random_state=42, max_iter=100)\n",
    "\n",
    "# fit the clustering to the the data\n",
    "kmeans.fit(data_skewed)\n",
    "\n",
    "# obtain the cluster labels\n",
    "kmeans_labels_skewed = kmeans.labels_\n",
    "\n",
    "# obtain the cluster means\n",
    "kmeans_means_skewed = kmeans.cluster_centers_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now visualize the results of the clustering as well as the obtained clusters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot petal length vs. petal width and corresponding classes\n",
    "scatter = ax.scatter(data_skewed[:,0], data_skewed[:,1], c=kmeans_labels_skewed.astype(np.float), cmap=plt.cm.Set1)\n",
    "\n",
    "# prepare data legend\n",
    "ax.legend(*scatter.legend_elements(), loc='upper left', title='Cluster')\n",
    "\n",
    "# plot cluster means\n",
    "ax.scatter(kmeans_means_skewed[:,0], kmeans_means_skewed[:,1], marker='x', c='black', s=100)\n",
    "\n",
    "# iterate over distinct cluster means\n",
    "for i, mean_skewed in enumerate(kmeans_means_skewed):\n",
    "    \n",
    "    # determine max cluster point distance\n",
    "    cluster_radi = cdist(data_skewed[:, 0:2][kmeans_labels_skewed==i], [mean_skewed]).max()\n",
    "    \n",
    "    # plot cluster size\n",
    "    ax.add_patch(plt.Circle(mean_skewed, cluster_radi, fc='darkgrey', edgecolor='slategrey', lw=1, alpha=0.1, zorder=1))\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[feature $x_1$]\", fontsize=14)\n",
    "ax.set_ylabel(\"[feature $x_2$]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Sample Skewed Distribution - k-Means Clustering Results', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the example above we observe one of the main disadvantages of the k-Means Clustering algorithm. k-Means Clustering has no way to account for oblong or elliptical clusters. Therefore, it fails to recognize the distinct clusters in very skewed distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. The Expectation Maximization (EM) Algorithm of Gaussian Mixtures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Expectation Maximization (EM)** algorithm technique is similar to the k-Means Clustering technique. However, instead of assigning examples to clusters to maximize the differences in means for continuous variables, the EM Clustering algorithm computes probabilities of cluster memberships based on one or more probability distributions. The goal of the clustering algorithm then is to maximize the overall probability or likelihood of the data, given the (final) clusters. Thereby, the objective of EM clustering is to estimate the means and standard deviations for each cluster so as to maximize the likelihood of the observed data (distribution). To achieve this objective the algorithm iteratively computes an (1) Expectation-Step and (2) Maximization-Step as described in the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step-1: The Expectation-Step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the **\"Expectation Step\"** we determine the expected values $E(z_{ij} | \\mu_{j}^{t}, \\sigma_{j}^{t})$ that denotes the probability of a given observation $x_i \\in \\mathcal{R}^d$ that it was drawn from the $j^{th}$ distribution (or cluster):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$E(z_{ij}) = \\frac{P(x_i|\\mu_{j}^{t}, \\sigma_{j}^{t})}{\\sum_{k'=1}^k P(x_i|\\mu_{k'}^{t},\\sigma_{k'}^{t})},$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "were:\n",
    "\n",
    "- $x_i \\in X$ denotes a particular data observation;\n",
    "- $t$ denotes the current iteration of the EM Clustering algorithm;\n",
    "- $\\mu_j$ and $\\sigma_j$ denotes the mean and (co-) variance of the $j^{th}$ distribution (or cluster)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the lecture you learned that the probability distribution of each cluster might be approximated by a Gaussian (Normal) probability distribution $\\mathcal{N}(\\mu, \\sigma)$. Hence the approach is usually referred to as the **EM-Algorithm of Gaussian Mixtures**. This simplification is justified by the application of the **\"law of large numbers\"** or **\"Central Limit Theorem\"** (you may want to have a look at the details of the theorem via the following link: https://en.wikipedia.org/wiki/Central_limit_theorem). In general, the probability density of a Gaussian \"Normal\" distribution, as defined by the formula below. It is parametrized its mean $\\mu$ and corresponding standard deviation $\\sigma$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ P(x|\\mu,\\sigma)  \\approx \\mathcal{N}(x|\\mu, \\sigma) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} e^{-\\frac{1}{2\\sigma^{2}}(x - \\mu)^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the **Central Limit Theorem** we can rewrite the formula of the Expectation Step as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$E(z_{ij}) = \\frac{P(x_i|\\mu_j^{t}, \\sigma_j^{t})}{\\sum_{k'=1}^k P(x_i|\\mu_{k'}^{t},\\sigma_{k'}^{t})} = \\frac{\\mathcal{N}(x_i|\\mu_j^{t}, \\sigma_j^{t})}{\\sum_{k'=1}^k \\mathcal{N}(x_i|\\mu_{k'}^{t},\\sigma_{k'}^{t})} = \\frac{ \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} e^{-\\frac{1}{2\\sigma_{j}^{t 2}}(x_i - \\mu_{j}^{t})^2} }{ \\sum_{k'=1}^{k} \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} e^{-\\frac{1}{2\\sigma_{k'}^{t 2}}(x_i - \\mu_{k'}^{t})^2} }$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distinct $E(z_{ij})$ values $z_{1j}$, $z_{2j}$, ..., $z_{nj}$ represent the probability distribution of the $j^{th}$ cluster that the data point is drawn from. Understanding the range of values the $z$ values can take is important. Each observation $x_i$ has k associated $E(z_{ij})$ values. In the **k-Means Clustering** algorithm each $z_{ij}$ can only take the value 0 or 1. This is why the k-Means Clustering algorithm is referred to as **\"hard\"** clustering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In contrast, the **\"EM-Clustering\"** algorithm is referred to as **\"soft\"** or **\"fuzzy\"** clustering. In EM-Clustering the distinct observations $x_i$ are considered to be drawn probabilistically from the distinct cluster distributions $j$. The corresponding $z_{ij}$ values can therefore be $z_{i1}$=0.85, $z_{i2}$=0.10 and $z_{i3}$=0.05, which represents a strong probability that the $x_{i}$ value originates from distribution (or cluster) 1 and a exhibit a smaller probability that it originates from distribution (or cluster) 2 or 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step-2: The Maximization-Step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the **\"Maximization-Step\"** we calculate update the parameters of each Gaussian \"Normal\" cluster distribution. Therefore, we derive for each $j^{th}$ distribution (or cluster) a new $\\mu_{j}$ and $\\sigma_{j}$ parameter value as given by:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\mu_{k}^{t+1} = \\frac {\\sum_{k'=1}^{k}E(z_{ij})x_i} {\\sum_{k'=1}^{k}E(z_{ij})}; \\sigma_{j}^{2, t+1} = \\frac {\\sum_{k'=1}^{k}E(z_{ij}) (x_i - \\mu_{k}^{t})^{2}} {\\sum_{k'=1}^{k}E(z_{ij})}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Expectation Maximization (EM) in a 2-Dimensional Feature Space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, to k-Means Clustering let's now define the parameters of the **EM-Clustering** algorithm. We will start by specifying the **number of clusters** $k$ we aim to detect in the iris dataset. We again hypothesize that our observations are drawn from an unknown distributions of three iris flower species (each distribution corresponding to a different mean $\\mu_1$, $\\mu_2$, and, $\\mu_3$ and corresponding standard deviation $\\sigma_1$, $\\sigma_2$, and, $\\sigma_3$). Therefore, we set the number of clusters to be detected to $k=3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_clusters = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to define a corresponding number of **initial 'means' $\\mu_{i}$** (the initial random cluster centers) that will be used as 'starting points' in the first iteration of the clustering process. In our case we will specify $k=3$ cluster means each of dimension 2, since we aim to retrieve 3 clusters based on the 2 features $x_1$ and $x_2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_means = np.array([[1.0, 3.0], [2.0, 6.0], [1.0, 7.0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will define a **maximum number of iterations** that we want to run the **EM-Clustering** algorithm. Please, note that the clustering terminates once there will be no further changes in the cluster assignments. However, it's good practice to define an upper bound of the iterations applied in the clustering (especially when analyzing datasets that exhibt a high-dimensional feature space):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_iterations = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we ready to initialize an instance of the **EM-Clustering** algorithm using Python's `sklearn` library of data science algorithms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "em = GaussianMixture(n_components=no_clusters, means_init=init_means, max_iter=max_iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the **EM-Clustering** to learn a model of the $x_1$ and $x_2$ features. Similar to the k-Means Clustering we will again make use of the `fit()` method provided by `sklearn` for each of its classifiers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "em.fit(data_skewed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have conducted the clustering, let's inspect the distinct cluster labels that have been assigned to the individual records of our artificial dataset. This can be achieved by calling the `predict()` function of the fitted model: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "em_labels_skewed = em.predict(data_skewed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, we want to inspect the coordinates of the cluster means (sometimes also referred to as \"centroids\") assigned by the algorithm. This can be achieved by calling the `means_`attribute of the fitted model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "em_means_skewed = em.means_ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now visually inspect the clustering results of the two features $x_1$ and $x_2$ in terms of the cluster assignments of each observation and cluster means learned by the **EM-Clustering** algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot x_1 vs. x_2 and corresponding cluster labels\n",
    "scatter = ax.scatter(data_skewed[:,0], data_skewed[:,1], c=em_labels_skewed.astype(np.float), cmap=plt.cm.Set1)\n",
    "\n",
    "# plot cluster means\n",
    "ax.scatter(em_means_skewed[:,0], em_means_skewed[:,1], marker='x', c='black', s=100)\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[feature $x_1$]\", fontsize=14)\n",
    "ax.set_ylabel(\"[feature $x_2$]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Sample Skewed Distribution - Expectation Maximization Results', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, we can observe that the EM Clustering technique nicely clusters the distinct distributions of the skewed distributions. We notice that the results are quite different from those computed by the k-Means clustering. As initially stated, the EM Clustering determines the **probability of cluster memberships** for each observation based on the initially specified $k$ probability distributions. Let's obtain those probabilites for each observation using the `predict_proba()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "em_probabilities_skewed = em.predict_proba(data_skewed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now visually inspect the clustering results of the two features $x_1$ and $x_2$ in terms of the learned probabilities of each observation corresponding to a particular cluster as determined by the **EM-Clustering** algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# init the plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "# add grid\n",
    "ax.grid(linestyle='dotted')\n",
    "\n",
    "# plot x_1 vs. x_2 and corresponding cluster labels\n",
    "scatter = ax.scatter(data_skewed[:,0], data_skewed[:,1], c=em_probabilities_skewed.astype(np.float), cmap=plt.cm.Set1)\n",
    "\n",
    "# plot cluster means\n",
    "ax.scatter(em_means_skewed[:,0], em_means_skewed[:,1], marker='x', c='black', s=100)\n",
    "\n",
    "# add axis legends\n",
    "ax.set_xlabel(\"[feature $x_1$]\", fontsize=14)\n",
    "ax.set_ylabel(\"[feature $x_2$]\", fontsize=14)\n",
    "\n",
    "# add plot title\n",
    "plt.title('Sample Skewed Distribution - Expectation Maximization Results', fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be observed that the EM-Clustering algorithm conducts a **\"soft\"** or **\"fuzzy\"** cluster assignment especially for observations that reside at the edges of the distinct clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lab Summary:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this fourth lab, a step by step introduction into the unsupervised **k-Means Clustering** and **EM-Clustering** machine learning algorithms was presented. The code and exercises presented in this lab may serve as a starting point for more complex and tailored programs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may want to execute the content of your lab outside of the Jupyter notebook environment, e.g. on a compute node or a server. The cell below converts the lab notebook into a standalone and executable python script. Pls. note that to convert the notebook, you need to install Python's **nbconvert** library and its extensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# installing the nbconvert library\n",
    "!pip3 install nbconvert\n",
    "!pip3 install jupyter_contrib_nbextensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now convert the Jupyter notebook into a plain Python script:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!jupyter nbconvert --to script cfds_lab_07.ipynb"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "1035.5999755859375px",
    "left": "38px",
    "top": "109.39990234375px",
    "width": "276.79998779296875px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
